%% Load the data - Select the File: 'ViscosityProfiles','TurningRate_QuasiSteady','TurningRate_Transient'
close all; clear all; clc;
Dir = {};
Dir{1} = '.\ViscosityProfiles\0pt5PEO';
Dir{2} = '.\ViscosityProfiles\0pt82PEO';
%
first = 4; last = 8; % Middle third of the channel
t = [0 10 20 30 40 50 60];
gradCalc_Start = 2; % 10 min
gradCalc_fin = 7; % 60 min
t_10x_5 = [5 15 25]; % time points to be filled in
x = []; Viscosity_profile = {};Grad_Eta_label = [];
for i = 1:length(Dir)
    load(Dir{i}); xlim_max = 1000; x(i,:) = xPos + (500 - median(xPos));
    Viscosity_profile{i} = cat(1,mu{1,:});
    
    FitViscSlope = [];FitVisc = [];
    for j = 1:length(mu)
        Mu_Exp = []; p = [];
        Mu_Exp = flip(mu{1,j}(first:last))
        p = polyfit(xPos(first:last),Mu_Exp,1)
        FitViscSlope(j,1) = p(1); FitViscSlope(j,2) = p(2);
        FitVisc = polyval(FitViscSlope(j,:),xPos(first:last)); % Create line of fitted data points
    end
    Grad_Eta_label(i,:) = mean(FitViscSlope(4:10,1)) % How the gradients are defined in Fig 1

    Grad_Eta(i,:) = FitViscSlope(gradCalc_Start:gradCalc_fin,1)
    % Find the center spatial viscosity (\eta_0)
    eta_0(i,:) = Viscosity_profile{i}(gradCalc_Start:gradCalc_fin,6);

    grad_eta_eta0(i,:) = Grad_Eta(i,:)./eta_0(i,:);
    p = [];
    p = polyfit(t(gradCalc_Start:gradCalc_fin),grad_eta_eta0(i,:),2);
    p_fit(i,:) = p(3)+p(2).*t(gradCalc_Start:gradCalc_fin)+p(1).*t(gradCalc_Start:gradCalc_fin).^2;

    p_10x_5(i,:) = p(3)+p(2).*t_10x_5+p(1).*t_10x_5.^2; % Calculate grad eta/eta_0 for 5, 15, 25 min
    grad_eta_eta0_5minInterval(i,:) = [p_10x_5(i,1) grad_eta_eta0(i,1) p_10x_5(i,2) grad_eta_eta0(i,2) p_10x_5(i,3)]
end
% Load the Data for short and WT
% Transient WT 5 cP:
load('.\TurningRate_Transient\Turning_Rates_WT_Transient_5(33-66).mat')
% Transient WT 9 cP (0.82% PEO):
load('.\TurningRate_Transient\Turning_Rates_WT_Transient_82(33-66).mat')
% Transient Short 5 cP (0.5% PEO):
load('.\TurningRate_Transient\Turning_Rates_ST_Transient_5(33-66).mat')
% Load steady-state short-flagella for different viscosity gradients
load('.\TurningRate_QuasiSteady\Turning_Rates_ST_(33-66)')

% Load and plot the steady-state wild-type for different viscosity gradients
close all;
load('.\TurningRate_QuasiSteady\Turning_Rates(33-66).mat')

close all;
MarkerSizeTransient = 3;
MarkerSizeSteady = 6;
MarkerBorderSteady = 1;
colors = {[0 0 0], [0.588 0.286 0.612], [0.929 0.694 0.125], [0.466 0.674 0.188], [0.85 0.325 0.0980], [0.301 0.745 0.933]};
colors_ST = {[0.301, 0.745, 0.933], [0.301, 0.745, 0.933], [0.301, 0.745, 0.933], [0.301, 0.745, 0.933], [0.301, 0.745, 0.933]};

grad_eta = [0.001, 0.79, 2.2, 3.4, 7.2]; grad_eta_std = [0.097827, 0.075021, 0.3408 ,0.96049, 1.9];

% Calculate grad eta/eta0 for 1.0% PEO:
GradEta_eta0 = grad_eta./center_vis;
PerPEO = [0 0.1 0.25 0.5 0.82];
p = polyfit(log(PerPEO(2:end)), log(GradEta_eta0(2:end)), 1);
m = p(1); b = exp(p(2)); x = linspace(0,1,100);
y_fit = b*x.^m; y_1PerPEO = b*1.^m;

GradEta_eta0_Short = [GradEta_eta0(1:4), y_1PerPEO];

Figure2d = figure('color',[1 1 1]); hold on; box on;
    % Plot steady-state short-flagella
    errorbarGradient_ST = [];
    for j = 1:length(GradEta_eta0_Short)
        colors_st = colors_ST{j};
        
        if j < length(GradEta_eta0_Short)
            y2 = (ones(100,1).*mean_omega_ST(j));
            x2 = [];
            temp_x = grad_eta(j)./center_vis(j);
            temp = sqrt((center_vis_std(j)./center_vis(j)).^2 + ones(1,1).*(grad_eta_std(j)/grad_eta(j)).^2);
            errorbarGradient_ST(j) = temp*temp_x;
            x2 = horzcat(x2, linspace(temp_x - temp*temp_x, temp_x + temp*temp_x,100)');
        else % Because 1.0% PEO gradient has no x-axis errorbars
           x2 = [];y2 = []; 
        end
        plot(x2,y2,'-','color',colors_st,'linew',MarkerBorderSteady)
        errorbar(GradEta_eta0_Short(j),mean_omega_ST(j),std_omega_ST(j),'^','color',colors_st,...
            'markerfacecolor','w', 'markersize', 8,'CapSize',0,'MarkerSize',MarkerSizeSteady,'LineWidth', MarkerBorderSteady)
    end
    
    % Plot fit to steady-state short-flagella
    [xData1, yData1] = prepareCurveData( GradEta_eta0_Short, mean_omega_ST );
    ft1 = fittype( {'x'}, 'independent', 'x', 'dependent', 'y', 'coefficients', {'a'} );
    [fitresult1, ~] = fit( xData1, yData1, ft1 );
    x_fit = linspace(0, 10, 100); y_fit1 = coeffvalues(fitresult1).*x_fit;
    p2 = plot(x_fit,y_fit1,'k--','linew',2);
    uistack(p2,'bottom');

    color_face = {[0.301, 0.745, 0.933],[0.466 0.674 0.188],[0.85 0.325 0.098]};
    color_edge = {[0.301, 0.745, 0.933],[0.466 0.674 0.188],[0.85 0.325 0.098]};
    % Plot transient data for short-flagella and wild-type
    f(1) = errorbar(grad_eta_eta0_5minInterval(1,:)*1000, mean_omega_ST_5,std_omega_ST_5, '^','color',color_edge{1},'MarkerSize',MarkerSizeTransient,'MarkerFaceColor',color_face{1},'CapSize',0);
    f(2) = errorbar(grad_eta_eta0_5minInterval(1,:)*1000, mean_omega_5,std_omega_5, 'o','color',color_edge{2},'MarkerSize',MarkerSizeTransient,'MarkerFaceColor',color_face{2},'CapSize',0);
    f(3) = errorbar(grad_eta_eta0_5minInterval(2,:)*1000, mean_omega_82,std_omega_82, 'o','color',color_edge{3},'MarkerSize',MarkerSizeTransient,'MarkerFaceColor',color_face{3},'CapSize',0);    
    
    % Plot steady-state wild-type
    errorbarGradient = [];
    for i = 1:length(grad_eta)
        y1 = []; x2 = [];
        x1 = ones(100,1)*grad_eta(i)./[center_vis(i)];
        y2 = (ones(100,1).*mean_omega(i));
        temp_x = grad_eta(i)./center_vis(i);  
        temp = sqrt((center_vis_std(i)./center_vis(i)).^2 + ones(1,1).*(grad_eta_std(i)/grad_eta(i)).^2);
        errorbarGradient(i) = temp*temp_x;
        y1 = horzcat(y1, linspace(mean_omega(i)-std_omega(i),mean_omega(i)+std_omega(i),100)'); % y-axis error boar
        x2 = horzcat(x2, linspace(temp_x - temp*temp_x, temp_x + temp*temp_x,100)'); % x-axis error boar
        % Error bars
        plot(x1,y1,'-','color',colors{i},'linew', MarkerBorderSteady)
        plot(x2,y2,'-','color',colors{i},'linew', MarkerBorderSteady)
        % Marker
        plot(grad_eta(i)./center_vis(i), mean_omega(i), 'o','markerfacecolor','w', 'markersize', MarkerSizeSteady, 'linew', MarkerBorderSteady, 'color', colors{i})
    
        grad_eta_eta(i) = (grad_eta(i))./center_vis(i);
        Tau(i) =  mean_omega(i); % Omega_visc
    end
    % Plot fit to wild-type
    [xData, yData] = prepareCurveData( grad_eta_eta, Tau );
    ft = fittype( {'x'}, 'independent', 'x', 'dependent', 'y', 'coefficients', {'a'} );
    [fitresult, ~] = fit( xData, yData, ft );
    y_fit = coeffvalues(fitresult).*x_fit;
    p1 = plot(x_fit,y_fit,'k','linew',2);
    uistack(p1,'bottom');
    hold off;

    xlim([0 4.7]); ylim([0 0.195]);
    xticks([0,1,2,3,4]);yticks([0,0.05,0.1,0.15]);
    xlabel('Normalized viscosity gradient \nabla\eta/\eta_0 (\mum^-^1\times10^-^3)');
    ylabel('Turning rate amplitude \omega_v_i_s_c (s^-^1)');
% Save Data
clc;
dir_excel = '.\Source data\MainText\';
tab = 4;
% Steady-state wild-type
writematrix(grad_eta_eta',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','A4')
writematrix(Tau',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','B4')
writematrix(errorbarGradient',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','C4')
writematrix(std_omega',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','D4')

% Steady-state short-flagella
writematrix(GradEta_eta0_Short',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','F4')
writematrix(mean_omega_ST',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','G4')
writematrix(errorbarGradient_ST',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','H4')
writematrix(std_omega_ST',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','I4')

% Transient wild-type: 3.4E-3 cP/um				
writematrix((grad_eta_eta0_5minInterval(1,:)*1000)',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','K4')
writematrix(mean_omega_5',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','L4')
writematrix(std_omega_5',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','M4')

% Transient wild-type: 7.2E-3 cP/um				
writematrix((grad_eta_eta0_5minInterval(2,:)*1000)',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','O4')
writematrix(mean_omega_82',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','P4')
writematrix(std_omega_82',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','Q4')

% Transient short-flagella: 3.4E-3 cP/um				
writematrix((grad_eta_eta0_5minInterval(1,:)*1000)',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','S4')
writematrix(mean_omega_ST_5',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','T4')
writematrix(std_omega_ST_5',[dir_excel,'Figure2.xlsx'],'Sheet',tab,'Range','U4')

alw = 0.5; %Border line width
mypapersize = [2 2]; %[width height]
set(gca,'LineWidth',alw);
set(gca,'XColor',[0 0 0]);set(gca,'YColor',[0 0 0]);
set(gcf, 'PaperUnits', 'inches', 'papersize', mypapersize, 'PaperPosition', [0 0 mypapersize]);

xticklabels({'','','','',''});
yticklabels({'','','',''});
ylabel('');xlabel('');title('')
ax = gca;outerpos = ax.OuterPosition;ti = ax.TightInset;
left = outerpos(1);bottom = outerpos(2);ax_width = outerpos(3);ax_height = outerpos(4);
ax.Position = [left+0.035 bottom+0.01 ax_width-0.04 ax_height-0.02];

% Save Figure
dir = '.\Figure2';
print(Figure2d , '-painters','-dpdf', '-r1200', [dir,'\Fig2d.pdf']); %r1200